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Abstract. - We consider a system in a non-equilibrium steady state by joining two semi-infinite 
Ising chains coupled to thermal reservoirs with different temperatures, T and T' . To compute the 
energy flux from the hot bath through our system into the cold bath, we exploit Glauber heat-bath 
dynamics to derive an exact equation for the two-spin correlations, which we solve for T' — oo 
and arbitrary T. We find that, in the T 1 — oo sector, the in-flux occurs only at the first spin. In 
the T < oo sector (sites x = 1,2, ...), the out-flux shows a non-trivial profile: F(x). Far from the 
junction of the two chains, F(x) decays as e~ x '^, where £ is twice the correlation length of the 
equilibrium Ising chain. As T — > 0, this decay crosses over to a power law (x~ 3 ) and resembles a 
"critical" system. Simulations affirm our analytic results. 



Introduction. — Equilibrium statistical mechanics, 
i ^ fo rmulated by Boltzmann, Gibbs, and others about a cen- 
tury ago, is so well established that today, it is part of the 
core material in typical physics programs. By contrast, 



statistical mechanics of systems far from equilibrium is so 



poorly understood that the Committee on CMMP 2010 
of the National Research Council recently recognized it 
lf~) as one of the six most fundamental and important chal- 
qq lenges in condensed matter and materials physics pQ . To 
f^be clear, we should distinguish two types of such sys- 
terns. On the one hand, we may encounter a statistical 
System evolving with a dynamics which obeys detailed 
;> balance. Though it will eventually settle into an equilib- 
J^~j rium state and can be described within the Boltzmann- 
Cibbs framework, its time- dependent behavior can be 
highly non-trivial. Good examples are slow relaxation or 
practically-stationary states, such as aging [2] or glassy 
phenomena J3J. On a more difficult level, we may face a 
system with a dynamics which violates detailed balance. 
Examples abound in, e.g., the life sciences, sociology, or 
economics, for which macroscopic variables are more ap- 
propriate. Typically, the rules for their evolution do not 
obey time-reversal. Even when such a system settles into 
a time-independent steady state, its stationary distribu- 
tion is unknown in general. Much of the intuition built 
on equilibrium statistical mechanics fails under these cir- 



cumstances. In an effort to develop insight into this class 
of non-equilibrium steady states (NESS), we take model 
systems with especially simple equilibrium properties and 
drive them out of equilibrium by mechanisms which are 
physically well motivated. 

A good example of a simple model driven far from equi- 
librium is the kinetic Ising model, coupled to two thermal 
reservoirs set at different temperatures [4 8. In previous 
studies, translational invariance is kept as much as possi- 
ble, e.g., every other spin being updated according to the 
same reservoir 06]. As a result, the effects of being driven 
out of equilibrium are extensive. In this article, we study 
another, perhaps more common, way of coupling a sys- 
tem to two baths, namely, one side being kept hot and the 
other, cold. Specifically, we consider an Ising chain (i.e., 
a one-dimensional lattice with sites x = — 1, 0, 1, ...) in 
which all spins with x < are coupled to one bath while 
the rest are coupled to another bath. Alternatively, we 
can view this system as joining together two Ising chains, 
each being in contact with its own thermal bath. The 
coupling of the two spins at the ends (x = 0, 1) allows the 
flow of energy from the hotter bath to the cooler one, as 
illustrated in fig. Q] 

Recently, a simpler version of this model has been in- 
vestigated Hj. On a ring of 2N + 1 spins (i.e., periodic 
chain), all but one is coupled to a standard bath at zero 
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Fig. 1: We consider two semi-infinite chains of Ising spins (red 
and blue) that only assume the two values a x — ±1, which 
we denote by "+" and "— " signs. The two boxes labeled T' 
and T, with T' > T, represent the two thermal reservoirs to 
which we couple the two Ising spin chains. The red (blue) 
arrows denote the energy flux flowing into (out of) the chain. 
Narrower arrows denote a smaller flux. 



temperature, while one spin is nipped randomly, mainly 
through a Poisson process. Energy is injected at the one 
spin, in that its bonds are randomly broken. Meanwhile, 
these broken bonds, or domain walls, wander into the rest 
of the ring and can annihilate each other when they meet. 
Since the main interest is how the details of the Poisson 
nipping control the fluctuations of injected energy, the fo- 
cus turned to the statistics of the location of the domain 
wall nearest the injection site. Thus, some information 
concerning the energy density profile is extracted. Since 
the domain wall statistics are the same as the statistics 
of bonds being broken or not along the chain, the domain 
wall information can be extracted from the correlation of 
nearest neighbor spins. Our work here differs in mainly 
two aspects. We find the entire two-spin correlation (not 
just the nearest neighbor pairs) and indicate how that re- 
sult provides us with correlations for arbitrary numbers 
of spins. Secondly, our approach allows for the study of 
the system set at any two temperatures, and we find ex- 
plicit results for the case where one part of the system is 
at infinite temperature and the other is at an arbitrary 
temperature. 

In the next section, the specifications of our model will 
be presented. The system is no longer translationally in- 
variant and the non-equilibrium effects are localized. The 
third section will be devoted to theoretical analysis and 
comparisons with simulation results. By restricting our- 
selves to "heat-bath" dynamics [TU|, we are able to solve 
for the correlation functions analytically. With these, the 
energy flux from the hotter bath to the colder one can be 
investigated using the approach in [5] . Since such an effect 
is localized (around x = 0), there is a non-trivial "pro- 
file" of the average energy flow between the baths and the 
chain. The details of this analysis are quite involved and 
will be published elsewhere [IT] . Here, we will only present 
the setup, some key steps, and the results. In particular, 
for the most extreme case: T x< o = oo and T x >o = 0, this 
profile vanishes for x < and decays as l/x 3 for x 3> 1. 
We end with a summary and outlook for future studies. 

The model. — In a kinetic Ising chain, spins are lo- 
cated on sites (labeled by an integer x) of a discrete lattice 
and assume only two values: o~ x = ±1. Thus, a configura- 



tion of this system is described by the set of spins {a x }. 
Each spin interacts with its two neighbors via the Hamil- 
tonian: H ({<J X }) = — JJ2x a xO~x+i, where J is a coupling 
constant. For a finite system, boundary conditions are 
needed, the simplest being periodic. 

While in contact with a thermal bath at temperature 
T, the kinetic Ising chain has a simple equilibrium dis- 
tribution: P cq ({&x}) = e _w / fcsT /Z. Its static proper- 
ties (e.g., the partition function Z) were obtained first 
by Ising [12 . Its dynamical properties, on the other 
hand, are much more complex since they depend sensi- 
tively on how the system is endowed with time depen- 
dence. The simplest is Glauber spin-flip "heat-bath" dy- 
namics [TU]. In a Monte Carlo simulation with parallel 
updates, each spin (a x ) is chosen and reset to be o~ x , with 
a probability that depends only on the average of its two 
neighboring spins: p x = \ + "fa x (cr x _i + a x +i) /4, with 
7 = tanh(2 J/ksT). Thus, with J, 7 > 0, the new spin 
will tend to be more aligned ((1 + j)/2) with its neighbors 
than anti-aligncd ((1 — 7)/2). The standard framework to 
deal with dynamics is based on the master equation for 
a time-dependent P ({o~ x },t). In this approach, we have 
P({a x },t + 1) = J2{a x }U x PxP({^x},t), and at large 
times, the system settles into the stationary P eq ({o~ x }). 
A favorite alternative update is random sequential, where 
only one randomly chosen spin is updated at each time 
step. We will use this method for the "two-temperature" 
Ising chain below. 

One way to drive this kinetic Ising chain out of equilib- 
rium is to couple it to more than one thermal bath. The 
simplest generalization is to have two baths, at T and 
X" > T. Even with this simplification, there is an infinite 
variety of ways to couple the spins to them. One is to 
update each spin the same way, randomly choosing T and 
T' with fixed probabilities [3] . On average, 7 remains ho- 
mogeneous and the intuitive picture of the system - same 
as in equilibrium but with an effective 7 - proves to be 
correct. One may alternatively couple the two baths to 
every other spin [5145l[T5] . so that we should write 7^ in 
the p x above. The master equation then takes the form 

P {{a x },t + 1)=^W (K} ; {a x }) P ({a x }, t) , (1) 



where W is given by 
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Here, -/V to t is the total number of spins in the system, so 
that N^l J2 X accounts for the possibility of any spin being 
chosen, with equal probability, once per time step. The 
riy^x factor insures that all other spins remain unchanged. 
For the model in [5], we have, say, 72,1 = tanh(2 J/fc^T) 
and 72n+i = tanh(2 J/ksT') with integer n. A promi- 
nent aspect of this model is the non-trivial energy flux 
through the system, flowing from the hotter bath to the 



p-2 



Energy flux near the junction of two Ising chains at different temperatures 



colder one. In particular, many of the properties of the 
energy transfered to the baths due to the updating of the 
spins can be computed exactly [B][B| . Entropy production 
(associated with the two baths) can be defined and also 
computed. Due to the translationally invariant nature of 
the alternating couplings, both the total flux and entropy 
production are extensive. 

In this study, we consider another way to couple the 
Ising chain to two baths, namely, 

j x < = tanh(2J/fc B T') = 7' (3) 
j x >! = tanh(2 J/k B T) = 7 (4) 

In other words, each of the two sectors of the chain is up- 
dated with a single temperature, with the left sector being 
hotter (T" > T). This is a much more common form of 
coupling a system to two baths, occurring, for example, 
in stovetop cooking. Of course, this model is still far from 
being realistic, since heat is typically transported, via dif- 
fusion, from the hotter side to the colder one. Here, due 
to the local nature of the spin-bath coupling, non-trivial 
energy transfer is expected to take place close to the junc- 
tion and the anomalies associated with non-equilibrium 
statistical mechanics should vanish for \x\ ^ 1. Now, 
translational invariance is broken and the analysis is more 
complex. In the next section, we provide a solution to 
the problem, in terms of the two-spin correlation function 
(a x o v ) . Simulation data will also be presented for the 
case T' = 00 (7' = 0), showing a consistent picture of the 
inhomogencous energy flux. 

Before ending, let us remark that we will, for simplicity, 
ignore the boundary conditions. A more rigorous way is 
to start with open boundaries, write modified expressions 
for the end spins, and ensure that they make no difference 
in the limit of large N to t ■ Alternatively, we can impose pe- 
riodic boundary conditions, deal with two junctions, and 
check that the effects of the two disentangle appropriately. 

Theoretical analysis and simulation results. 

For W's which violate detailed balance, it is generally 
impossible to find an explicit expression for the stationary 
distribution P* {{<r x }), i.e. the solution to P* = WP*. 
A formal solution exists [14] , but it is so cumbersome that 
typically, useful information cannot be extracted. An al- 
ternative approach is to seek the correlations functions in 
the steady state 

(a Xl . . . a Xk ) =^2a Xl ... a Xk P* ({a x }) . (5) 

Inserting this definition into the equation for P*, we face 
the sum J2{a x } °"xi ■ • • Vx k W ({a x }; {a x }) which, in gen- 
eral, involves k' > k spins (er's). Writing an equation for 
the fc'-spin correlation, we encounter more spins. The re- 
sult is the BBGKY hierarchy [T5] involving correlations 
of arbitrarily many spins. However, W is linear in its 
arguments in our case (a salient feature of "heat-bath" 




Fig. 2: The domains R , R++, and in which we solve 

eq. [10] for (a x a y ). We have the BC {a x a x ) = 1 along the red 
line and match the solutions along the green lines. 

dynamics), so that coupled equations do not proliferate. 
In particular, we find 

a " W = Y + fro+O + (-^tot - 1) ffa • (6) 

{**} 

We now replace the dummy a a by a x and average both 
sides of eq. [6] The A to t variable cancels and we obtain 

= ~f x (er x _i + a x+ i) - 2 (a x ) (7) 

(with slight modification for the spins at the end of the 
chain). Unless j x = 1, the only solution to this equation is 
(a x ) = 0. This result is hardly surprising, especially since 
an ordinary Ising chain displays no spontaneous magneti- 
zation for all positive temperatures. 

Turning to the two-point function, we first note 

(a x a x ) = 1 (8) 

and 

(a x a y ) = (a y a x ) (9) 
so that we can focus on, say, x < y only. There, we obtain 

= 7z ((Vx + lVy) + (<T x -l<Ty)) 

+ 7 V ((<t x <t v+ i) + (<J x o- y -i)) - 4 (cr x a y ) . (10) 

It is easy to verify that, for "f x = ^y y = 7, a familiar result 
emerges, i.e., (c x (J y ) cq = w'^ - ^, with w = tanh ( J/ksT). 

Before we solve equation eq. ([TTT[) , let us remark that all 
correlations with an odd number of spins vanish. Further, 
with our dynamics, correlations of any even number of 
spins can be expressed in terms of products of {(J x a y ) 
[T1H3JHS1IIZI • I n tms sense, we have a complete picture of 
this non-equilibrium system. 

Returning to our problem, we see it is a generalized, 
discrete Hclmholtz equation and so, boundary conditions 
(BCs) arc crucial. One BC is eq. [8] and the other is 
(a x a y ) — > for \x — y\ — > 00. In other words, we have 
Dirichlet BCs for the half space x > y. In our case, we 
have only two 7's. Thus, this domain can be partitioned 
into three regions, R , and (corresponding to 
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x, y being positive or negative, as shown in fig. [2]), in each 
of which the operator is homogeneous. Thus, the solution 
can be attacked by standard, though tedious, methods. 
We will only indicate some key points for the general prob- 
lem and provide a few more details below for the special 
case of T" = oo, i.e., 7' = 0. 

Note that the operator in cq. (fit))) is isotropic in R++ 

and R , but anisotropic in Since the operators are 

homogeneous (discrete) Laplacians within each region, we 
can make use of the eigenfunctions e lkx + l py with eigenval- 
ues 2 (j x cos fc — 1) + 2 (7 y cosp — 1). In the regions R 

and R++, appropriate combinations must be chosen to 
satisfy eqs. ([HIO- The last BC is less straightforward and 
we will rely on the following approximation. We choose 
functions which vanish on the x = N line in and the 

y = —N line in R first, and then let TV — >• 00 at the end. 

Such an approach can also be exploited in R^ Finally, 

we must match the functions on the common boundaries 
(x = and y = lines, i.e. the green lines in fig. E]) - 

Let us consider an easier case, 7' = 0, so we have only 
one parameter left: 7. Eq. (|10p immediately provides 



(o~ x a y ) — for y < x < 



(11) 



which is quite reasonable for R where spins are nipped 

completely randomly. Next, in , this equation simpli- 
fies to a single one, (a x a y+ i) + (a x a y -i) = 4 (o- x a y ) /j, for 
all y < (and x > y). It is satisfied by (cr x a y ) = e^ x , with 
(real) fi = ±cosh _1 (2/7). To be careful, let us first focus 
on y < 0. From eq. (fTTj) above, we have (aoo~y<o) = 0, 
which restricts us to (a x a y ) = Asmh/ix. Imposing the 
BC for x — > 00, we arrive at 



{<J x o-y<o) = V x > y . 



(12) 



On the line y = 0, we have a different BC, since (coOo) 
1. The other BC then picks out 



(cr x a ) = uj x 



V x > 



(13) 



where 



1 



(14) 



7 \Z(t, 

In other words, this particular correlation is of the same 
form as that for the Ising chain in equilibrium, with an ef- 
fective temperature given by 7/2 instead of 7. Intuitively, 
this result is appealing, since one spin (00) is coupled to 
7' = and the other (a y ) is coupled to 7. 

Finally, we turn our attention to R++. Here, eq. (jTUJ) 
is exactly the same as for the equilibrium Ising chain. In- 
deed, the only difference between the two problems is the 
BC on the y = line, where we have eq. ([T3")) instead of 
(c x <7o) cq = lo x . Thus, we can simplify our problem by 



considering the difference 

S (x,y) = (< 



(15) 



which vanishes on two of the three boundaries of R++. A 
further advantage is S is directly related to the quantity 



of interest - the flux, F (x), i.e., the average rate of energy 
loss when o~ x >a is updated. Specifically, from [B], one can 
show that 

F{x) = - [jS (x + 1, x - 1) - S (x + 1, x) - S (x, x - 1)] 

(16) 

in units of 4J. In other words, once we find S, the energy 
flux profile can be computed. 

Our problem now reduces to solving the Dirichlet prob- 
lem 

1 



VS = ; 
in R ++ with BC's 

S (x, 0) = Cj x - oj' 1 



V 



1 



7 



(17) 



S(x,x) = S(N,y)=0 (18) 



Here, A x , y is the discrete Laplacian and N will be taken 
to 00 at the end. A standard route is to use the Dirichlet 
Green's function G which satisfies G = on the bound- 
aries and T>G(x, y; £, 77) = S x .^S ViT] . The result is 

G (x v £ n) = ^ V V Uk *( x >V y > Uk * &lH (19) 
^,y,^,n) N2 2^ 2^ 2- 7 (cosfc + cosp) ' 1 > 

where the eigenfunctions of our operator (T> and BC's) are 

Uk,p (£j v) — sinfc^ sinp?7 — sinp£ sinfcjy (20) 

with (k,p) = ir{m,n)/N. Exploiting the identity S = 
^[ST^G — GVS] and a discrete divergence theorem, we 
find 

N-l 

S(x, y) = 2\2 [G(x, y; f , 1) (ufi - ^)] . (21) 

Applying 4/7 = cj + cj -1 and 2/7 = cj + cj _1 , we arrive at 
an explicit solution in R++ : 



N 2 



sinfcsinp (cosfc — cosp) 
2 — 7(cos k + cosp) 

[7 (cos k + cosp) — 3] 
(2 — 7 cosfc) (2 — 7 cosp) (1 — 7 cosfc) (1 — 7 cosp) 

(22) 

In the N — > 00 limit, p, fc become continuous in [0, n] 
and convenient variables are 6 = (fc + p)/2 and <f> = 
(fc — p) /2. Inserting the result into the expression for F, 
the flux profile can be written as a double integral 

3 /-7T [ir-4> 

ir( x ) = J_ / / d9d(f>A{9,(j))sm29x, (23) 
^ Jo Jo 

where x = 1, 2, . . . and 

sin 8 sin 2 <p [sin 2 (j> — sin 2 0] [7 cos <fi — cos 9] 



A l 



1 — 7 cos 9 cos < 
(-1)« 



4^ 2a 2 - 47a cos 8 cos cj) + 7 2 (cos 29 + cos 2<fi) 
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We now consider the large x asymptotics of F. Defer- 
ring a detailed presentation to elsewhere jllj . we provide 
only some results here. We find that eq. [231 f° r integers 
x > 2 reduces exactly to 



(1 - tV) (l - - 7V) 



2.r 



F(x) = 5— rrr I — ; 7 dtl, 

(25) 

with a known correction for x = 1 To get the cross- 
over behavior as 7 — > 1 (T — > 0), we can arrange an 
asymptotic expansion of eq. [25] (for large x) as a sum over 
modified Bessel functions 
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a = 0.997 ±0.001 
b= -5,92 ±0.01 



: 0.997 ± 0.1103 
6 = -12.60 ±0.04 



= W 2 * £] B„ C" n e C ^n+i (0 . 



(26) 



n=0 



where £ = 2yl — 7 2 x and B„ are explicitly known coeffi- 
cients. 

Intuitively we expect F to decay exponentially for T > 
0, where the equilibrium correlation length, £ = — l/lnw, 
is finite. Careful analysis confirms this picture. In 
particular, eq. ([2"rJ]) leads to 
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Fig. 3: The average (over 10 4 independent simulation runs) of 
the change in energy AE and the variance (AE) 2 — AE at 
spin x = 3 as a function of the number of times t that the spin 
was updated per run. Here, 7 = 0.9 (with T' = 00). 
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Interestingly, we observe that the decay length is £/2, 
though we have no good heuristic explanation. In the 
extreme case of T = 0, these complex expressions simplify 
and a power law appears. The result is 



lim F(x) 

7— >i 



1 



illustrated in fig. [21 both the average and the variance in- 
crease with t linearly (i.e., a = I). The intercept, b, of the 
former provides us with data to compare with F(x). Sim- 
(I + O^- 1 )}. (27) ilarly, the intercept of the latter allows us to estimate the 
error bars in fig. [4] Fig. [4] shows the excellent agreement 
with theoretical predictions and simulations. 

Finally, we remark that we did not present results for 
the in- flux of energy, which occurs only at x = 0, since 
there are no correlations involving a x< Q. It is clear that 
its value must be just T, x> qF (x). 



2ttx 3 



(28) 



Though reminiscent of the crossover to critical exponents, 
we also have no simple argument in favor of the specific 
value —3. Of course, the total steady flux is finite, so that 
T, x> qF (x) must converge. But this constraint only limits 
the power to be less than —2. 

As further confirmation of these analytic results, we 
performed Monte-Carlo simulations of our system. For 
T 1 = 00, it is unnecessary to include any spin at x < 0. 
Instead, we study a chain with N + 1 spins, the first (ctq) 
being coupled to the infinite temperature bath. At each 
step, a random spin is chosen and assigned a new value ac- 
cording to the probabilities given above, except <To , which 
is assigned ±1 randomly. Using N w 30, we evolved the 
system for approximately 10 6 time steps (i.e. updates per 
spin). When a spin flips during an update, we record the 
change in the energy AE, which can take only the values 
and ±1 (in units of 4J). In this manner, we have a time se- 
ries AE (x, t) for each site. It is clear that, over a run, this 
quantity essentially performs a random walk. For systems 
in equilibrium, such walks will be unbiased. For our case, 
we expect biased walks. To verify this scenario, we carried 
out 10 4 independent runs and computed both the average, 

9 r 1 2 

AE (x, t), and the variance [AE(x,t)\ - \AE(x,t)\ . As 



Summary and Outlook. — In this article, we con- 
sider joining together two Ising chains end-to-end, each 
coupled to a uniform (but different) temperature reser- 
voir. As a result, when the system settles into a steady 
state, there is a constant flow of energy through the sys- 
tem, from the hotter to the colder bath. Far from the 
junction, each sector of the chain is expected to behave 
as one in equilibrium, so that the effects of this energy 
flow arc localized. Such a through-flux is strictly a non- 
equilibrium phenomenon and the stationary distribution 
is not simply given by some spatially dependent, effective 
Boltzmann factor. By choosing heat-bath spin-flip dynam- 
ics, we are able to compute all the correlation functions. 
In particular, we provide key points for the calculation 
and the results for the two-spin correlation, (a x a y ). This 
is not simply a function of \x — y\, since our system is not 
translationally invariant. From these functions, we can 
compute F(x), the steady flux of energy out of the sys- 
tem, as a result of updating the spin at x. When one 
reservoir is set at T" = 00, we find explicitly that F (x) 
decays into the sector with T < 00 exponentially, with 
a length half of the correlation length of its equilibrium 
counterpart. As T — > 0, there is no such counterpart, as 
the decay crosses over to a power law: 1/x 3 . Simulations 
of this model verify these results. 
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Fig. 4: The dots are the energy flux F(x) at spin x computed 
via Monte Carlo simulations (with error bars smaller than the 
dots for 7 > 0.9). The solid lines are the first two orders of the 
expansion in eq. [26] 



Before closing, let us consider possible avenues for fu- 
ture studies, as well as related Ising-likc systems far from 
equilibrium. The most immediate case of interest is fi- 
nite T' , analytic results for which should be within reach. 
There, we expect exponentials for F (x) in both sectors. It 
is not clear if the decay lengths continue to be as simple as 
what we found. In particular, F must vanish as T" — > T. 
Whether this limit appears as an overall amplitude which 
vanishes, analytically with (T" — T) 2 or more singularly, 
would be interesting. Beyond the Ising chain, the most 
obvious generalization is the same model in higher dimen- 
sions, joining two systems side by side, say. Since the ordi- 
nary Ising model at equilibrium displays phase transitions 
in higher dimensions, more interesting questions arise. Ex- 
amples include: Does F decay with singular powers when 
T is set at T c , the critical temperature (with T" T c )? If 
so, are there new exponents (such as the dynamic expo- 
nent z in model A j!8))? or new combinations of existing 
exponents (such as z = 4 — rj in model B [TS])? Set- 
ting T < T c , we could expect a non-trivial magnetization 
profile (a x ). How does it decay with xl Continuing along 
these lines, are there new singularities associated with this 
profile when T" is lowered to T c l 

Another generalization of joining two Ising systems is 
also immediate: using Kawasaki spin-exchange dynamics 
instead. While no new static properties are present 
in an equilibrium Ising model with this type of dynamics, 
much more surprising behavior for the two-temperature 
case have been observed, for both 1- and 2-d systems [5D] . 
Beyond these, there are limitless varieties, most of which 
will undoubtedly provide further insight into how we may 
proceed in attempting to formulate an overarching princi- 
ple for non-equilibrium statistical mechanics. 
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